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ABSTRACT 

We study predictions for dark matter phase-space structure near the Sun based on high- 
resolution simulations of six galaxy halos taken from the Aquarius Project. The local DM den- 
sity distribution is predicted to be remarkably smooth; the density at the Sun differs from the 
mean over a best-fit ellipsoidal equidensity contour by less than 15% at the 99.9% confidence 
level. The local velocity distribution is also very smooth, but it differs systematically from a 
(multivariate) Gaussian distribution. This is not due to the presence of individual clumps or 
streams, but to broad features in the velocity modulus and energy distributions that are stable 
both in space and time and reflect the detailed assembly history of each halo. These features 
have a significant impact on the signals predicted for WIMP and axion searches. For example, 
WIMP recoil rates can deviate by ~ 10% from those expected from the best-fit multivariate 
Gaussian models. The axion spectra in our simulations typically peak at lower frequencies 
than in the case of multivariate Gaussian velocity distributions. Also in this case, the spec- 
tra show significant imprints of the formation of the halo. This implies that once direct DM 
detection has become routine, features in the detector signal will allow us to study the dark 
matter assembly history of the Milky Way. A new field, "dark matter astronomy", will then 
emerge. 
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1 INTRODUCTION 

In the 75 years since Zwicky ( 1933) first pointed out the need 
for substantial amounts of unseen material in the Coma cluster, 
the case for a gravitationally dominant component of non-baryonic 
dark matt er has become overwhelmingly strong. It seemed a long 
shot when lPeeblesI dl982h first suggested that the dark matter might 
be an entirely new, weakly interacting, neutral particle with very 
low thermal velocities in the early universe, but such Cold Dark 
Matter (CDM) is now generally regarded as the most plausible 
and consistent identification for the dark matter. Particle physics 
has suggested many possible CDM particles beyond the standard 

model. Two promising candidates are WIMPs (weakly interact- 
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ing massive particles, see Lee & Weinberg 1977; Gunn et al. 1978; 

Ellis et al. 1984) and axions jPeccei & Ouinnl l 1977b a: W einberg! 
1 97 8l: I Wilczekl 1 1 97 8h ■ Among the WIMPs, the lightest supersym- 
metric particle, the neutralino, is currently favoured as the most 
likely CDM particle, and the case will be enormously strength- 
ened if the LHC confirms supersymmetry. However, ultimate con- 



firmation of the CDM paradigm can only come through the di- 
rect or indirect detection of the CDM particles themselves. Neu- 
tralinos, for example, are their own antiparticles and can annihi- 
late to produce 7-rays and other particles. One goal of the recently 
laun ched Fermi Gamma-ray space telescope is to dete ct this radia- 
tion dGehrels & Michelsonlll999l:ISpringel et alj2008h . 

Direct detection experiments, on the other hand, search for 
the interaction of CDM particles with laboratory apparatus. For 
WIMPs, detection is based on nuclear recoil events in mas- 
sive, cryogenically co oled bolometers in underground laboratories 
dJungman et ai] [l996); for axions, resonant microwave cavities in 
strong magnetic fi elds exploit the axion-photon conversion pro- 
cess isikivid 1985b . Despite intensive searche s, the only experimen t 
which has so far reported a signal is DAMA (Bernabe i et alj|2007h 
which has clear evidence for an annual modulation of their event 
rate of the kind expected from the Earth's motion around the Sun. 
The interpretation of this result is controversial, since it appears 
to require dark matter properties which are in conflict with up- 
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per limits established by other experimen ts (see lSavage et alj2004l ; 
iGondolo & Gelminil200l ; lGelminil200^ . for a discussion and pos- 
sible solutions). Regardless of this, recent improvements in detector 
technology may enable a detection of "standard model" WIMPS or 
axions within a few years. 

Event rates in all direct detection experiments are determined 
by the local DM phase-space distribution at the Earth's position. 
The relevant scales are those of the apparatus and so are extremely 
small from an astronomical point of view. As a result, interpret- 
ing null results as excluding specific regions of candidate param- 
eter space must rely on (strong) assumptions about the fine-scale 
structure of phase-space in the inner Galaxy. In most analyses the 
dark matter has been assumed to be smoothly and spherically dis- 
tributed about the Gala ctic Centre with an isotropic Maxwellian ve- 
locity distribution (e.g. Freese et al. 1988) or a mult i variate Gaus- 
sian distribution (e.g. lUllio & Kamionkowskill200ll ; lGreenll200ll ; 
iHelmi et alj|2002h . The theoretical justification for these assump- 
tions is weak, and when numerical simulations of halo formation 
reached sufficiently high resolution, it became clear that the phase- 
space of CDM halos contains considerable substructure, both grav- 
itationally bound subhalos and unbound streams. As numerical res- 
olution has improved it has become possible to see structure closer 
and closer to the centre, and this has led some investigators to sug- 
gest that the CDM distribution near the Sun could, in fact, be almost 
fractal, with large density variations over short length-scales (e.g. 
iKamionkowski & Koushiappasll2008l) . This would have substantial 
consequences for the ability of direct detection experiments to con- 
strain particle properties. 

Until very recently, simulation studies were unable to resolve 
any subs tructure in regions as close to the Galactic Centre as the 
Sun (see lMoore et all 2001; Hel mi et al . 2002l l2003l, for example). 
This prevented realistic evaluation of the likelihood that massive 
streams, clumps or holes in the dark matter distribution could af- 
fect event rates in Earth-bound detectors and so weaken the par- 
ticle plry_si£s_j^ojiclujiiOTS 

(see lSavage et alj|2006l : IKa mionkows ki & Koush iappas 20081 f° r 
recent discussions). As we shall show in this paper, a new age has 
dawned. As part of its Aquarius Project dSpringel et al .l[2008h the 
Virgo Consortium has carried out a suite of ultra-high resolution 
simulations of a series of Milky Way-sized CDM halos. Simula- 
tions of ind ividual Milky Way ha los o f similar scale have been car- 
ried out bv lDiemand et alj d2008h and lStadel et"ail ( 2008). Here we 
use the Aquarius simulations to provide the first reliable character- 
isations of the local dark matter phase-space distribution and of the 
detector signals which should be anticipated in WIMP and axion 
searches. 



2 THE NUMERICAL SIMULATIONS 

The cosmological parameters for the Aquarius simulation set are 
fi m = 0.25, Ov = 0.75,(T8 = 0.9, n s = 1 and Ho = 
100 h km s _1 Mpc -1 with h — 0.73, where all quantities have 
their standard definitions. These parameters are consistent with cur- 
rent cosmological constraints within their uncertainties, in partic- 
ular, with the parame ters inferred from the WMAP 1-year an d 
5-year data analyses dSpergel et all 120031 ; iKomatsu et"al] 120081) . 
Milky Way-like halos were selected for resimulation from a par- 
ent cosmological simulation which used 900 3 particles to follow 
the dark matter distribution in a 100/i -1 Mpc periodic box. Se- 
lection was based primarily on halo mass (~ 10 12 M Q ) but also 
required that there should be no close and massive neighbour at 
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Figure 1. Top panel: Density probability distribution function (DPDF) for 
all resimulations of halo Aq-A measured within a thick ellipsoidal shell 
between equidensity surfaces with major axes of 6 and 12 kpc. The lo- 
cal dark matter density at the position of each particle, estimated using an 
SPH smoothing technique, is divided by the density of the best-fit, ellip- 
soidally stratified, power-law model. The DPDF gives the distribution of 
the local density in units of that predicted by the smooth model at random 
points within the ellipsoidal shell. At these radii only resolution levels 1 
and 2 are sufficient to follow substructure. As a result, the characteristic 
power-law tail due to subhalos is not visible at lower resolution. The fluc- 
tuation distribution of the smooth component is dominated by noise in our 
64-particle SPH density estimates. The density distribution measured for 
a uniform (Poisson) particle distribution is indicated by the black dashed 
line. Bottom panel: As above, but for all level-2 halos after rescaling to 
Vmax = 208.49 km/s. In all cases the core of the DPDF is dominated by 
measurement noise and the fraction of points in the power law tail due to 
subhalos is very small. The chance that the Sun lies within a subhalo is 
~ 10~ 4 . With high probability the local density is close to the mean value 
averaged over the Sun's ellipsoidal shell. 



z = 0. The Aquarius Project resimulated six such halos at a series 
of higher resolutions. The naming convention uses the tags Aq-A 
through Aq-F to refer to these six halos. An additional suffix 1 to 
5 denotes the resolution level. Aq-A-1 is the highest resolution cal- 
culation, with a particle mass of 1.712 x 10 3 M Q and a virial mass 
of 1.839 x 10 12 M© it has more than a billion particles within the 
virial radius R200 which we define as the radius containing a mean 
density 200 times the critical value. The Plummer equivalent soft- 
ening length of this run is 20.5 pc. Level-2 simulations are available 
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for all six halos with about 200 million particles within ifeoo- Fur- 
ther details of the ha los and their characteristics can be found in 
ISpringeletaP l l2008h . 

In the following analysis we will often compare the six level-2 
resolution halos, Aq-A-2 to Aq-F-2. To facilitate this comparison, 
we scale the halos in mass and radius by the constant required to 
give each a maximum circular velocity of V ma x = 208.49 km/s, 
the value for Aq-A-2. We will also sometimes refer to a coordi- 
nate system that is aligned with the principal axes of the inner halo, 
and which labels particles by an ellipsoidal radius r e ii defined as 
the semi-major axis length of the ellipsoidal equidensity surface on 
which the particle sits. We determine the orientation and shape of 
these ellipsoids as follows. For each halo we begin by diagonal- 
ising the moment of inertia tensor of the dark matter within the 
spherical shell 6 kpc < r < 12 kpc (after scaling to a com- 
mon Vmax). This gives us a first estimate of the orientation and 
shape of the best fitting ellipsoid. We then reselect particles with 
6 kpc < r c n < 12 kpc, recalculate the moment of inertia tensor 
and repeat until convergence. The resulting ellipsoids have minor- 
to-major axis ratios which vary from 0.39 for Aq-B-2 to 0.59 for 
Aq-D-2. The radius restriction reflects our desire to probe the dark 
matter distribution near the Sun. 



3 SPATIAL DISTRIBUTIONS 

The density of DM particles at the Earth determines the flux of 
DM particles passing through laboratory detectors. It is important, 
therefore, to determine not only the mean value of the DM density 
8 kpc from the Galactic Centre, but also the fluctuations around this 
mean which may result from small-scale structure. 

We estimate the local DM distribution at each point in our 
simulations using an SPH smoothing kernel adapted to the 64 
nearest neighbours. We then fit a power law to the resulting dis- 
tribution of lnp against lnr c n over the ellipsoidal radius range 
6 kpc < r e ii < 12 kpc. This defines a smooth model density 
field Pmodei(feii). We then construct a density probability distribu- 
tion function (DPDF) as the histogram of p/p mo dei for all particles 
in 6 kpc < r e u < 12 kpc, where each is weighted by p _1 so that 
the resulting distribution refers to random points within our ellip- 
soidal shell rather than to random mass elements. We normalise the 
resulting DPDFs to have unit integral. They then provide a prob- 
ability distribution for the local dark matter density at a random 
point in units of that predicted by the best fitting smooth ellipsoidal 
model. 

In Fig. Q] we show the DPDFs measured in this way for all 
resimulations of Aq-A (top panel) and for all level-2 halos after 
scaling to a common V max (bottom panel). Two distinct compo- 
nents are evident in both plots. One is smoothly and log-normally 
distributed around p = p mo dci, the other is a power-law tail to high 
densities which contains less than 10~ 4 of all points. The power- 
law tail is not present in the lower resolution halos (Aq-A-3, Aq- 
A-4, Aq-A-5) because they are unable to resolve subhalos in these 
inner regions. However, Aq-A-2 and Aq-A-1 give quite similar re- 
sults, suggesting that resolution level 2 is sufficient to get a reason- 
able estimate of the overall level of the tail. A comparison of the six 
level 2 simulations then demonstrates that this tail has similar shape 
in different halos, but a normalisation which can vary by a factor 
of several. In none of our halos does the fraction of the distribu- 
tion in this tail rise above 5 x 10~ 5 . Furthermore, the arguments of 
Springel et al (2008) suggest that the total mass fraction in the in- 
ner halo (and thus also the total volume fraction) in subhalos below 
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Figure 2. Top four panels: Velocity distributions in a 2 kpc box at the Solar 
Circle for halo Aq-A-1. ui, V2 and V3 are the velocity components parallel 
to the major, intermediate and minor axes of the velocity ellipsoid; v is the 
modulus of the velocity vector. Red lines show the histograms measured 
directly from the simulation, while black dashed lines show a multivari- 
ate Gaussian model fit to the individual component distributions. Residuals 
from this model are shown in the upper part of each panel. The major axis 
velocity distribution is clearly platykurtic, whereas the other two distribu- 
tions are leptokurtic. All three are very smooth, showing no evidence for 
spikes due to individual streams. In contrast, the distribution of the velocity 
modulus, shown in the upper left panel, shows broad bumps and dips with 
amplitudes of up to ten percent of the distribution maximum. Lower panel: 
Velocity modulus distributions for all 2 kpc boxes centred between 7 and 
9 kpc from the centre of Aq-A- 1 . At each velocity a thick red line gives the 
median of all the measured distributions, while a dashed black line gives 
the median of all the fitted multivariate Gaussians. The dark and light blue 
contours enclose 68% and 95% of all the measured distributions at each ve- 
locity. The bumps seen in the distribution for a single box are clearly present 
with similar amplitude in all boxes, and so also in the median curve. The 
bin size is 5 km/s in all plots. 



4 Vogelsberger et al. 




v [km s '] v [km s 1 ] 

5f 1 1 1 ] 5r 1 1 




v [km s '] v [km s 1 ] 

Figure 3. Distributions of the velocity modulus in four well separated 2 kpc 
boxes about 8 kpc from the centre of Aq-A. Results are shown for each 
region from each of the three highest resolution simulations. Error bars are 
based on Poisson statistics. The different resolutions agree within their error 
bars, and show the same bumps in all four boxes. For the purpose of this 
plot, we have chosen a larger bin for our histograms, 10 km/s as compared 
to 5 km/s in our other velocity plots. For this bin size the statistical noise in 
Aq-A- 1 is barely visible. 



the Aq-A-1 resolution limit is at most about equal to that above this 
limit. Hence, the chance that the Sun resides in a bound subhalo of 
any mass is of order 10 -4 . 

The striking similarity of the smooth log-normal component 
in all the distributions of Fig.[T]has nothing to do with actual den- 
sity variations in the smooth dark matter distribution. It is, in fact, 
simply a reflection of the noise in our local density estimates. We 
demonstrate this by setting up a uniform Poisson point distribu- 
tion within a periodic box and then using an SPH smoothing kernel 
adapted to the 64 nearest neighbours to associate a local density 
with each particle in exactly the same way as for our halo simula- 
tions. We can then construct a DPDF for these estimates (relative to 
their mean) in exactly the same way as before. The result is shown 
in the top panel of Fig.Q]as a dashed black line. It is an almost per- 
fect fit to the smooth component in the simulations, and it would fit 
the other halos equally well if plotted in the lower panel. 

The fit is not perfect, however, and it is possible to disentan- 
gle the true scatter in density about the smooth model from the 
estimation noise. The latter is expected to be asymptotically log- 
normal for large neighbour numbers, and Fig. [T] shows that it is 
very close to log-normal for our chosen parameters. If we assume 
that the scatter in intrinsic density about the smooth model is also 
approximately log-normal, we can estimate its scatter as the square 
root of the difference between the variance of the simulation scat- 
ter and that of the noise: symbolically, a lnt v = \J o-'£ hs — o-'^ oisc . 
Indeed, it turns out that the variance in ln(p/p mo doi) which we 
measure for our simulated halos (excluding the power-law tail) is 
consistently higher than that which we find for our uniform Poisson 
distribution. Furthermore, tests show that the differences are stable 
if we change the number of neighbours used in the SPH estima- 
tor to 32 or 128, even though this changes the noise variance by 



factors of two. This procedure give the following estimates for rms 
intrinsic scatter around the smooth model density field in our six 
level-2 halos, Aq-A-2 to Aq-F-2: 2.2%, 4.4%, 3.7%, 2.1%, 4.9% 
and 4.0% respectively. The very large particle number in the radial 
range we analyse results in a standard error on these number which 
is well below 0.01 for all halos. Thus, we can say with better than 
99.9% confidence that the DM density at the Sun's position dif- 
fers by less than 15% from the average over the ellipsoidal shell 
on which the Sun sits. This small scatter implies that the density 
field in the inner halo is remarkably well described by a smooth, 
ellipsoidal, power-law model. 

We conclude that the local density distribution of dark matter 
should be very smooth. Bound clumps are very unlikely to have 
any effect on direct detection experiments. The main reason for 
this is the short dynamical time at the solar radius (about 1% of 
the Hubble time). This results in very efficient mixing of unbound 
material and the stripping of all initially bound objects to a small 
fraction of the maximu m mass they may have had in the past (see 
Vogelsberg er et alj2008l , for a discussion of these processes). Note 
that the actual density of DM in the Solar neighbourhood and the 
shape of the equidensity surfaces of the Milky Way's dark mat- 
ter distribution will depend on how the gravitational effects of 
the baryonic components have modified structure during the sys- 
tem's formation. Unfortunately, the shape of the inne r DM halo of 
the Milky Way is poorly constrained observationally jHelmill2004 
lLaw et alf2 005). The dissipative contraction of the visible compo- 
nents probably increased the density of the dark matte r component, 
and may also have made it m ore axisymmetric (e.g. lGnedin et al.l 
|2004 iKazantzidis et alj[2004l) but these processes are unlikely to 
affect the level of small-scale structure. The very smooth behaviour 
we find in our pure dark matter halos should apply also to the more 
complex real Milky Way. 



4 VELOCITY DISTRIBUTIONS 

The velocity distribution of DM particles near the Sun is also an 
important factor influencing the signal expected in direct detec- 
tion experiments. As mentioned in the Introduction, most previ- 
ous work has assumed this distribution to be smooth, and either 
Maxwellian or multivariate Gaussian. Very different distributions 
are possible in principle. For example, if the local density distribu- 
tion is a superposition of a relatively small number of DM streams, 
the local velocity distribution would be effectively dis crete with all 
particles in a given stream sharing the same v elocity dSikivie et al.l 
1 19951 : IStiff etalfeoOll : IStiff & Widrowll2003l) . Clearly, it is impor- 
tant to understand whether such a distribution is indeed expected, 
and whether a significant fraction of the local mass density could 
be part of any individual stream. 

We address this issue by dividing the inner regions of each of 
our halos into cubic boxes 2 kpc on a side, and focusing on those 
boxes centred between 7 and 9 kpc from halo centre. In Aq-A-1, 
each 2 kpc box contains 10 4 to 10 5 particles, while in the level- 
2 halos they contain an order of magnitude fewer. For every box 
we calculate a velocity dispersion tensor and study the distribu- 
tion of the velocity components along its principal axes. In almost 
all boxes these axes are closely aligned with those the ellipsoidal 
equidensity contours discussed in the last section. We also study 
the distribution of the modulus of the velocity vector within each 
box. The upper four panels of Fig. [2] show these distributions of 
a typical 2 kpc box at the solar circle in Aq-A-1 (solid red lines). 
Here and in the following plots we normalise distributions to have 
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Figure 4. Velocity modulus distributions in exactly the same format as the bottom panel of Fig. [2] but for all six of our halos at level-2 resolution. All 
distributions are smooth. Only in Aq-B-2 do we see a strong spike which is due to a single box which has 60% of its mass (though a small fraction of 
its volume) in a single subhalo. No other box in any of the distributions has a subhalo contributing more than 1.5% of the mass. All distributions show 
characteristic broad bumps which are present in all boxes in a given halo, and so in its median distribution. These bumps are in different places in different 
halos. 



unit integral. The black dashed lines in each panel show a multi- 
variate Gaussian distribution with the same mean and dispersion 
along each of the principal axes. The difference between the two 
distributions in each panel is plotted separately just above it, This 
particular box is quite typical, in that we almost always find the 
velocity distribution to be significantly anisotropic, with a major 



axis velocity distribution which is platykurtic, and distributions of 
the other two components which are leptokurtic. Thus the veloc- 
ity distribution differs significantly from Maxwellian, or even from 
a multivariate Gaussian. The individual velocity components have 
very smooth distributions with no sign of spikes due to individual 
streams. This also is a feature which is common to almost all our 
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Figure 5. Left panel: Median velocity modulus distributions for all six level-2 halos repeated from Fig. [4] The black dashed line is the mean of these 
distributions. Middle panel: Deviations of the velocity modulus distribution of each of the six halos from the sample mean. The amplitude of the various 
bumps is similar in different halos and over the whole velocity range. It reaches more than 10% of the amplitude of the mean distribution. Right panel: Relative 
deviations of the individual velocity modulus distributions from their sample mean. Typical relative deviations are about 30%, but they can exceed 50% at 
higher velocities. 



2 kpc boxes. It is thus surprising that the distribution of the velocity 
modulus shows clear features in the form of bumps and dips with 
amplitudes of several tens of percent. 

To see how these features vary with position, we overlaid the 
distributions of the velocity modulus for all 2 kpc boxes centred 
between 7 and 9 kpc from the centre of Aq-A- 1 (bottom panel of 
Fig.[2](. We superpose both the directly measured distributions and 
the predictions from the best-fit multivariate Gaussians. At each 
velocity, the solid red line show the median value of all the directly 
measured distributions, while the dashed black line is the median of 
all the multivariate Gaussian fits. The dark and light regions enclose 
68% and 95% of all the individual measured distributions at each 
velocity. 

It is interesting to note that the bumps in the velocity distri- 
bution occur at approximately the same velocity in all boxes. This 
suggests that they do not reflect local structures, but rather some 
global property of the inner halo. In Fig.[3]we show velocity modu- 
lus distribution for four different boxes in Aq-A at the three highest 
resolutions (levels 1, 2 and 3). The error bars are based on Poisson 
statistics in each velocity bin. Clearly the same bumps are present 
in all boxes and at all resolutions. Thus, they are a consequence of 
real dynamical structure that converges with increasing numerical 
resolution. 

In Fig. [4] we make similar plots of the velocity modulus dis- 
tribution for all level-2 halos. These distributions are quite smooth. 
The sharp peak in Aq-B-2 is due to a single 2 kpc box where 60% 
of the mass is contained in a single subhalo. No other box in this 
or any other halo has more than 1.5% of its mass in a single sub- 
halo. The great majority of boxes contain no resolved subhalo at 
all. Although the details of the median distributions vary between 
halos, they share some common features. The low velocity region 
is more strongly populated in all cases than predicted by the mul- 
tivariate Gaussian model. In all cases, the peak of the distribution 
is depressed relative to the multivariate Gaussian. At moderately 
high velocities there is typically an excess. Finally, and perhaps 
most importantly, all the distributions show bumps and dips of the 
kind discussed above. These features appear in different places in 
different halos, but they appear at similar places for all boxes in a 



given halo. The left panel of Fig. [5] superposes the median velocity 
modulus distributions of all level-2 halos and plots their mean as a 
black dashed line. The middle panel shows the deviations of the in- 
dividual halos from this mean. The amplitudes of the deviations are 
similar in different halos and at low and high velocities. In percent- 
age terms the deviations are largest at high velocity reaching values 
of 50% or more, as can be seen from the right panel of Fig. [5] 

The bumps in the velocity distribution are too broad to be 
explained by single streams. Furthermore, single streams are not 
massive enough to account for these features. This is shown more 
clearly in Fig. [6] where we illustrate some streams in velocity space 
for a 2 kpc box in halo Aq-A-1. Different colours here indicate 
particles that belonged to different FoF groups at redshift 4.2. For 
clarity we only show streams from groups that contribute at least 
10 particles to this volume (0.025% of the total number of particles 
present at this location). There are 27 such objects. If we consider 
all FoF groups that contribute more than 2 particles to the volume 
shown in Fig. [6] we find that a given FoF group contributes streams 
that are typically only populated by 2 particles (0.005% of the 
total mass in the box). This implies that most of the groups con- 
tribute several streams of very low density. The most prominent 
streams have ~ 40 particles, i.e. ~ 0.1% of the mass in this vol- 
ume. This clearly shows that streams are expected to be neither 
dense nor massive in the Solar vicinity. 

The most prominent streams typically occupy the tail of the 
velocity distribution in these local boxes. The excess of parti- 
cles moving with similar velocities and high kinetic energies can 
be measured using a velocity correlation function, as shown by 
lHeimietal.U2002h . 



5 ENERGY DISTRIBUTIONS 

We have seen that the distributions of individual velocity compo- 
nents in localised regions of space are very smooth, whereas the 
velocity modulus distribution shows clear bumps. Taken together 
with the fact that these bumps occur at similar velocities in regions 
on opposite sides of the halo centre, this indicates that we must be 
seeing features in the energy distribution of dark matter particles. 
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Figure 6. Streams in velocity space for a 2 kpc box ~ 8 kpc from the cen- 
tre of Aq-A- 1 . Different colours stand for particles associated to different 
FoF groups at redshift 4.2. Only groups contributing more than ten particles 
are shown. The box contains 27 such objects and has in total 41143 parti- 
cles (shown as small black points) of which 1796 come from these groups. 
Clearly, particles originating from the same group cluster in velocity space 
and build streams; often many streams per group. 



Figure 7. Mean phase-space density distribution as a function of energy for 
Aq-A for particles in a spherical shell between 6 and 12 kpc and for all 
five resolution levels. Especially at high-binding energies the convergence 
is very good. Features in the distribution function are visible at all resolu- 
tions for energies below 2.7 V^ ax , despite the fact that the mass resolution 
differs by more than a factor of 1800 between Aq-A-1 and Aq-A-5. The 
less bound parts show more variation from resolution to resolution but still 
agree well between Aq-A- 1 and Aq-A-2. 



To investigate this further, we estimate the mean phase-space 
density as a function of energy in each of our halos using the prop- 
erties of the particles at radii between 6 and 12 kpc. Clearly our 
halos are not perfectly in equilibrium and they are far from spher- 
ical. Thus their phase-space densities will only approximately be 
describable as functions of the integrals of motion, and they will de- 
pend significantly on integrals other than the energy. Nevertheless 
we can estimate a mean phase-space density as a function of energy 
by taking the total mass of particles with 6 kpc < r < 12 kpc and 
energies in some small interval and dividing it by the total phase- 
space volume corresponding to this radius and energy range, e.g. 

HF) - dM 1 m 
where f(E) is the energy-dependent mean phase-space density and 



j d 3 x s/2 (E- 



*(*)), 



(2) 



9(E) = 4t 

V,B>*(x) 

is the available phase-space volume in the configuration-space vol- 
ume V. The differential energy distribution is easily calculated by 
binning the energies of all particles between 6 and 12 kpc. The 
phase-space volume can be calculated by solving for the gravita- 
tional potential at the position of all simulation particles and then 
using these as a Monte-Carlo sampling of configuration space in the 
relevant integrals. Taking the ratio then yields the desired estimate 
of /(B). 

In Fig.[7Jwe show f(E) measured in this way for all our sim- 
ulations of Aq-A. We express the energy in units of Vm ax and we 
take the zero-point of the gravitational potential to be its average 
value on a sphere of radius 8 kpc. As a result the measured en- 
ergy distribution extends to slightly negative values. Note how well 
the distribution converges at the more strongly bound energies. At 
higher energies the convergence between the level- 1 and 2 reso- 
lutions is still very good. This demonstrates that we can robustly 
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Figure 8. Evolution of the mean phase-space density distribution of Aq-A- 
2 over four time intervals of about 300 Myr. Below 2.4 V^ ax the phase- 
space distribution function is time-independent, implying that the system 
has reached coarse-grained equilibrium. The small bumps at these energies 
are therefore well-mixed features in action space. The variability of the fea- 
tures in the weakly bound part of the distribution shows that they are due 
to individual streams and therefore change on the timescale of an orbital 
period. Note that the phase-space density at these energies is almost three 
orders of magnitude below that of the most bound particles. 



T 




Figure 9. Scaled phase-space distribution functions for all level-2 halos. 
In addition to scaling according to V m ax we have also corrected for a 
zero-point offset in the potential energy between different halos. The black 
dashed line shows the average distribution function based for our halo sam- 
ple. At high binding energies the scatter between average and individual 
halo distribution functions is quite small, showing that this part of the distri- 
bution function is near-universal. At low-binding energies large amplitude 
features are visible in all halos. These features differ from halo to halo and 
are related to recent events in their formation histories. 



measure the mean phase-space density distribution. Furthermore, 
we see clear wiggles that reproduce quite precisely between the 
different resolutions. 

Fig. [8] shows similarly estimated mean phase-space density 
distributions for Aq-A-2 at five different times separated by about 
300 Myr. This is longer than typical orbital periods in the region we 
are studying. Despite this, the wiggles at energies below 2.4 V^ ax 
are present over the complete redshift range shown. This demon- 
strates that these features are well-mixed, and that the phase-space 
distribution function has reached a coarse-grained equilibrium. In 
contrast, the variability of the wiggles in the part of the distribution 
corresponding to weakly bound particles (where the orbital periods 
are much larger) shows that these must be due to individual streams 
or to superpositions of small numbers of streams, which have not 
yet phase-mixed away. 

To estimate what these phase-space distribution functions 
should look like for a "smooth" system, we average the functions 
found in our six individual level-2 halos. In Fig. [9] we superpose 
these six functions and their mean (/) (the black dashed line). The 
similarity of the different distribution functions at high binding en- 
ergies suggests a near-universal shape for f(E). At lower binding 
energies, individual halos deviate quite strongly from (/). This can 
be seen more clearly in Fig.llOlwhere we plot log(//{/)), the dec- 
imal logarithm of the ratio of the phase-space density of an indi- 
vidual halo to the mean. The lower axis is orbital energy in units of 
Knax> while the upper axis is the corresponding dark matter parti- 
cle velocity at the Solar Circle. In this plot one can clearly see the 
wiggles, which are located at different energies for different halos. 
For Vs k P c < 350 km/s the distribution functions for all halos sat- 
isfy 0.7 < //(/) < 1.4. For low binding energies (velocities of 
600 km/s or more at the Solar Circle) this ratio can exceed a factor 
of ten. 

These features in the phase-space density distribution must be 
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Figure 10. Deviations of the individual phase-space density distributions 
from the mean over our sample of level-2 halos. We focus here on the more 
bound part. The lower x-axis shows the orbital energy while the upper one 
shows the corresponding velocity 8 kpc distance from halo centre. The am- 
plitude of features increases for Vg k pc > 350 km/s. At even lower binding 
energies, E > 3 V m ax deviations can reach an order of magnitude, see 
Fig.! 

related to events in the formation of each halo. To demonstrate this 
explicitly, we have computed f(E) separately for particles which 
were accreted onto two of our halos (i.e. first entered the main 
progenitor FoF group) at different epochs. The upper left panel of 
Fig. II II shows that Aq-A-2 had a very "quiet" merger history. Ma- 
terial accreted at different times is arranged in a very orderly way 
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Figure 11. Top row: Contributions to the present-day phase-space density distribution from particles accreted at different epochs (indicated by different 
colours). The top left panel shows the build-up of the distribution function for halo Aq-A-2. This halo has a quiescent formation history with no recent 
mergers. The top right panel is a similar plot for Aq-F-2, which underwent a major merger between z = 0.75 and z = 0.68. Bottom row: These plots isolate 
the contribution of a single, massive FoF group to the z = phase-space density distributions. For Aq-A-2 this group was identified at z = 6.85, for Aq-F-2 at 
z = 0.75. In both cases it is clear that material from the group is responsible for some of the features seen in the present-day phase-space density distribution. 



in energy space. All the most strongly bound particles were ac- 
creted before redshift 5, and material accreted at successively later 
times forms a series of "shells" in energy space. The most weakly 
bound wiggles are due entirely to the most recently accreted ma- 
terial, and progressively more bound bumps can be identified with 
material accreted at earlier and earlier times. In contrast, the top 
right panel shows that Aq-F-2 had a very "active" merger history, 
with a major merger between z = 0.75 and z = 0.68. The corre- 
spondence between binding energy and epoch of accretion is much 
less regular than for Aq-A-2, and much of the most bound material 
actually comes from the object which fell in between z — 0.75 
and z = 0.68. It is also striking that many of the wiggles in this 
object are present in material that accreted at quite different times, 
suggesting that they may be non-steady coherent oscillations rather 
than stable structures in energy space. Nevertheless, in both ha- 
los one can identify features in the phase-space density distribu- 
tion with particles accreted at certain epochs, and in both halos the 
most weakly bound particles were added only very recently. Note 
that the phase-space density of this material is very low, so it con- 
tributes negligibly to the overall local dark matter density. In the 
bottom panels of Fig.[TT]we show the f(E) distributions of parti- 
cles which were associated with a single, massive FoF-group which 



was identified at z = 6.85 in the case of Aq-A-2 and at z — 0.75 
in the case of Aq-F-2. The wiggles in the strongly bound part of 
Aq-A-2 are clearly due to this early merger event, while the later 
merger in Aq-F-2 is responsible for most of the material accreted 
in 0.56 < z < 0.81 and for most of the strong features in the 
phase-space density distribution. 

We conclude that these features in the energy distribution 
should open the window to "dark matter astronomy" once exper- 
iments reach the sensitivity needed for routine detection of DM 
particles. We will then be able to explore the formation history of 
the Milky Way using the DM energy distribution. 



6 DETECTOR SIGNALS 

We will now use the spatial and velocity distributions explored 
above to calculate expected detector signals. The main question 
here is how the non-Gaussian features of the velocity distribution 
influence these signals. Our results show that features due to sub- 
halos or massive streams are expected to be unimportant. On the 
other hand, deviations of the velocity distributions from a perfect 
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Figure 12. Recoil spectra ratio for the three highest resolution simulations of Aq-A (left) and the level-2 (right) simulations of the other halos. For these plots 
we averaged the recoil rate over a year for every box and then calculated the median recoil rate ratio (-R)/(-R g auss) of the rates for the simulation and for 
the best-fit multivariate Gaussian distribution. The x-axis is directly proportional to the energy. In all level-2 halos the expected recoil spectrum based on a 
multivariate Gaussian can be wrong by about 10% depending on the energy. Furthermore the behaviour of the deviations seems quite similar. This is due to 
the fact that the velocity distributions differ in a characteristic way from a multivariate Gaussian. The deviations in the recoil spectra are typically highest at 
high energies. 



Gaussians in terms of general shape, bumps and dips can have an 
impact on detector signals. 

There are currently more than 20 direct detection experiments 
searching for Galactic DM, most of them focusing on WIMPs. For 
these, the detection scheme is based on nuclear recoil with the de- 
tector mate rial. The differential WIMP elastic scattering rate can be 
written as I Jung man et alj[l99r3) : 

R = n Po T(E,t), (3) 

where 1Z encapsulates the particle physics parameters (mass and 
cross-section of the WIMP; form factor and mass of target nucleus), 
po is the local dark matter density that we assume to be constant 
based on the results of section 3, and 

oo 

T(E,t)= fdv^, (4) 

^min 

where f v is the WIMP speed distribution in the rest frame of the 
detector integrated over the angular distribution. v m in here is the 
detector-dependent minimum WIMP speed that can cause a recoil 
of energy E: 

where m x is the WIMP mass and tua the atomic mass of the target 
nucleus. To get detector independent results we set 1Z — 1 in the 
following^. 

The recoil rate shows a annual modulation over the year 
jDrukieretaT]|l986h . To take this into account we add the Earth's 
motion to the local box velocities to transform Galactic rest frame 
velocities into the detector frame. We model the motion of the 

1 This also implies that we assume the form factor to be constant. Any 
other form factor will change the shape of the recoil spectrum. Since we are 
not interested in the exact shape of the spectrum, but in deviations expected 
due to different velocity distributions, we neglect form factor effects in the 
following. 



Earth according to lLewin & Smith! dl996h and lBinnev & Merrifleldl 
( 1998). Let ve — u r + us + ue be the velocity of the Earth rel- 
ative to the Galactic rest frame decomposed into Galactic rotation 
u r , the Sun's peculiar motion us and the Earth's velocity relative to 
the Sun ue- In Galactic coordinates these velocities can be written 
as: 

u r = (0, 222.2, 0) km/s, 
us = (10.0,5.2,7.2) km/s, 
SE,i = mb(A) cos(/3i) sin(A - A;), 
Mb (A) = (he) [1 — esin(A — A,:)] (6) 

where i — R,cj>,z, A is the ecliptic longitude (Ao = (13 ± 
1)°),(ub) = 29.79 km/s is the mean velocity of the Earth around 
the Sun, and the ellipticity of the Earth orbit e = 0.016722. The 
u r value is based on a combination of a large number of indepen- 
dent d eterminations of the circular velocity bv lKerr & Lvnden-Belll 
(1986). We note that this value has a standard deviation of 20 km/s. 
For the constant /3 and A angles we take: 

(J3r,P<I>,0x) = (-5.5303°, 59.575°, 29.812°) 
(A R ,A ,A Z ) = (266.141°, -13.3485°, 179.3212°) (7) 

The ecliptic longitude can be written as 

X(t) = L{t) + 1.915° sin g(t) + 0.020° sin 2g(t), 
Lit) = 280.460° +0.9856474% 

git) = 357.528° + 0.9856003% (8) 

where t is the fractional day number relative to noon (UT) on 31 
December 1999 (J2000.0). We refer to a day number relative to 31 
December 2008 in our plots. In what follows we will assume that 
the indirection is always aligned with the major axis of the princi- 
pal axis frame of the velocity ellipsoid in each box and the (j>- and 
z-directions, with the intermediate and short axes. This is needed 
to add the Earth's motion to the box velocities, and to transform the 
velocity vectors in each box to the detector frame. 

Clearly the deviations of the velocity distribution from a per- 
fect multivariate Gaussian found in the previous sections will also 
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Figure 13. Top panels: Annual modulation for all 2 kpc boxes with halocentric distance between 7 and 9 kpc in halo Aq-A-1 assuming v mln = 300 km/s. 
The left plot shows how the dimensionless recoil rate (R(t) — (R))/(R) changes over the year. The right plot shows the corresponding modulation parameter 
space defined by the peak day (a>axis) and maximum amplitude (R, max — (R))/ (R) (2/ _ax i s )- Bottom panels: Modulation parameters for the local 2 kpc 
boxes of all level-2 resolution halos. There is no clear trend visible in the day of maximum behaviour over the halo sample. On the other hand, the median 
amplitude in all boxes is higher than expected based on the Gaussian sample for v mln = 300 km/s. The line and contour scheme is the same as in Fig. [5] 



alter the recoil spectrum, because the velocity integral T(E, i) ef- 
fectively measures the 1/ii-weighted area under the velocity curve. 
As in the previous sections we compare the results obtained di- 
rectly from the simulations to the expectation for a best-fit multi- 
variate Gaussian distribution. In Fig. [12] we plot recoil spectra ra- 
tios for the three highest resolution simulations of Aq-A (left) and 
the level-2 (right) simulations of the other halos. For these plots 
we averaged the recoil rate over a year for individual boxes. The 
rates are calculated using the simulation velocity distribution ((R)) 
and the best-fit Gaussian model for each box ((i? gauss )). The plots 
show the median of the ratios (R) / (R gsiuss ) over all boxes. Since 
we assume that the density po is constant in a given box, it drops 
out when calculating the ratios. The a;-axis measures the energy in 
dimensionless j3 = v/c values. For a given detector this can easily 
be converted to keV, assuming the masses m x and tua are given in 
GeV/c 2 : 



E : 



2m 2 x m A j ,j lnl , 



c /T x 10° keV 



(9) 



(m x + mA) 2 

Fig.[T2]clearly shows that in all level-2 halos the expected re 



coil spectrum based on a multivariate Gaussian model can differ 
by up to 10% from the directly predicted simulation result. Fur- 
thermore, the behaviour of the deviations seems to be similar in all 
cases, especially at low energies, where we already found that the 
phase-space density is nearly universal. The similarity in the devia- 
tions between the different halos is due to the fact that the velocity 
distributions all differ in a characteristic way from the Gaussian dis- 
tributions as shown in section 4. The deviations in the recoil spectra 
are typically highest at high energies. 

The 10% deviations in the recoil spectra are larger than the 
typical deviations expected due to the annual modulation. There- 
fore these deviations from the Gaussian model can also influence 
the annual modulation signal. In Fig. [T3] (top row) we plot the di- 
mensionless recoil rate (R(t) — (R) ) / (R) of all local 2 kpc boxes at 
~ 8 kpc from the centre of Aq-A-1 (left),where (R) is the annual 
mean recoil. We have assumed f m i n = 300 km/s for all plots in 
this figure. The maximum can clearly be seen around the month of 
June. The plot on the right in Fig.[l3]shows the modulation parame- 
ter space that we define by the day of maximum amplitude (x-axis) 
and the maximum modulation amplitude of the recoil rate over the 
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Figure 14. Top panel: Median recoil rate amplitude for all 2 kpc boxes 
with halocentric distance between 7 and 9 kpc for all level-2 halos. The 
plot shows the difference between the relative maximum modulation am- 
plitude observed in the simulation and that expected for the best-fit multi- 
variate Gaussian distribution. Bottom panel: Median day of maximum am- 
plitude for the same halos (solid red) compared to their Gaussian predic- 
tions (dashed black). The day of maximum amplitude is the same for all 
boxes and is well reproduced in the Gaussian model. The phase-reversal 
can clearly be seen. 



year defined as (i? max — (R))/{R) (y-axis). The bottom row of 
Fig.[]~3]shows the maximum amplitude (left) and day of maximum 
(right) for all level-2 halos (solid red) and the corresponding best-fit 
multivariate Gaussian model (dashed black). 

Comparing the Gaussian median values to the box median val- 
ues one can see that the day of maximum amplitude does not devi- 
ate significantly from that predicted for a multivariate Gaussian; in 
particular there is no clear trend visible over the halo sample. On 
the other hand, the median amplitude in all halos is slightly higher 
than expected based on the Gaussian sample for « m in = 300 km/s. 

The amplitude differences for various n m in values are shown 
in Fig. [14] Here we calculated the maximum amplitude and day of 
maximum for different n m in values for all level-2 halos. The ampli- 
tude plot (top) shows the difference between the maximum relative 
modulation amplitude observed in the simulation and that expected 
for the best-fit multivariate Gaussian model. The maximum ampli- 
tude Umin-dependence is similar for the six halos. Since only the 
velocity distribution enters into the recoil calculation, this similar- 
ity is due to the fact that deviations of the halo velocity distribution 



from the Gaussian model are also quite similar for all six halos. 
The bottom plot of Fig.[74lshows the day of maximum amplitude is 
well predicted by the multivariate Gaussian for all halos. The sharp 
transition in th e day of maximum is due to the well-known phase- 
reversal effect (Prim ack et al.ll 19881) . We checked that the subhalo 
dominated box in Aq-B-2, where by chance about 60% of the box 
mass is in a single subhalo, leads to a very different modulation sig- 
nal. The day of maximum in that case shifts about 100 days from 
the Gaussian distribution. We note that although the subhalo mass 
fraction in this particular box is high, the subhalo volume fraction 
is tiny, so even within this box, almost all observers would see the 
smooth regular signal. 

Although most of the direct experiments currently search 
for WIMPs the axion provides another promising candidate for 
CDM. It arises from the Peccei-Quinn solution to the strong CP- 
problem. One axion detection scheme is based on using the axion- 
electromagnetic coupling to induce resonant conversions of axions 
to photons in the microwave frequency range. Galactic axions have 
non-relativistic velocities (J3 = v/c ~ 10 -3 ) and the axion-to- 
photon conversion process conserves energy, so that the frequency 
of converted photons can be written as: 



Au a = 241.8 



lAteV/c 2 



l + i/3 2 ] MHz (10) 



where m a is the axion mass that lies between 10 6 eV/c^ and 



10 



-3 



eV/c 2 . 5/xeV axions would therefore convert into v„ = 



1200 MHz photons with an upward spread of A ~= 2 kHz due 
to their kinetic energy. An advantage of axion detection compared 
to WIMP searches is the fact that it is directly sensitive to the en- 
ergy rather than to the integral over the velocity distribution. The 
power P, developed in the axion search cavity due to resonant 
axion-photon conversion can be written as dSikividl 1983b : 



P = V pa cavity) 



(ID 



where V encapsulates the experimental properties (cavity vol- 
ume, magnetic field, quality factor) and particle physics properties 
(model dependent coupling parameter, axion mass). The only as- 
trophysical input is the local density p a (y ca vit y ) of axions with en- 
ergies corresponding to the cavity frequency. For simplicity we set 
V = 1. We can produce axion spectra from our simulations by tak- 
ing a local volume element (a box) and computing the distribution 
of kinetic energies K of the particles found in this location. The 
number of particles with a given K is then directly proportional to 
p a at this frequency, and so to the power in the frequency bin. 

To make the results independent of axion mass and other ex- 
perimental properties we present histograms of (3 2 normalised to 
one. For a given axion mass m a (in ^teV/c 2 ) the x-axis must be 
transformed according to x — * 241.8 m a (1 + l/2i) to get the 
corresponding frequencies in MHz. 

A long-running a xion search experiment is ADMX at LLNL 
dHagmann et al.| [l996). It has channels at medium (MR) and high 
resolution (HR). The latter has a frequency resolution of about 
0.02 Hz. For u° = 500 MHz and an axion velocity of v = 
200 km/s this translates into a velocity width of only 0.018 km/fl 
Our numerical resolution prevents us from predicting the behaviour 
on such small scales. For wider bin searches and especially for the 
medium resolution (MR) channel (125 Hz corresponding to a typi- 
cal velocity spread of about 100 km/s) we can, however make reli- 
able predictions by binning particles with respect to /3 2 . 



For non-relativistic motion we can write dv = (c 2 /v) (dv/v®). 
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Figure 15. Axion spectra of level-2 halos for all 2 kpc boxes with halocentric distance between 7 and 9 kpc. Rescaling the x-axis according to x — > 
241.8 m a (1 + l/2x) for an axion mass m a in /ieW yields the x-axis in MHz. The y-axis is proportional to the power P developed in the detector cavity. 
Therefore the panels show the predicted frequency spectra expected for an axion search experiment like ADMX. These spectra can be reasonably described 
by a multivariate Gaussian but significant differences remain. The maximum in the power is at lower frequencies in the simulation than in the Gaussian model. 
The bumps already found in the velocity and energy distribution are clearly visible in these spectra. In all halos the power at low frequencies is higher than 
expected from the Gaussian model. The line and contour scheme are the same as in Fig. [2] 

In Fig. [T5] we show axion spectra for all level-2 halofl In a broad sense, the spectra obtained from our simulations look similar 



We neglect the effects of the Earth's motion when constructing the spec- cally leads to a shift of about 100 Hz due to annual modulation and a daily 

tra since here our focus is on the general spectral shape. This motion typi- shift of about 1 Hz due to the Earth's rotation JPuffv et al.l2005t) . 
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to those of multivariate Gaussian models. However there are a num- 
ber of differences. For example, the peak power is shifted to lower 
frequencies. The Gaussian distribution is also a poor description of 
the spectrum at low frequencies. In all halos the power at low and 
high frequencies is higher then expected from a multivariate Gaus- 
sian model. This effect is quite small for high frequencies but very 
significant for low frequencies. The higher power at low frequen- 
cies can be understood from the velocity distributions in Fig. [4] 
In Aq-B-2 the subhalo dominated box that was seen in Fig. [4] is 
clearly visible as a peak in the power spectrum at high frequency. 
The bumps in the velocity distribution also result in quite signifi- 
cant features in the axion spectra that might be visible in the MR 
channel given enough signal-to-noise. 



7 CONCLUSION 

We have characterised the local phase-space distribution of dark 
matter using the recently published ultra-high resolution simula- 
tions of the Aquarius Project. Our study provides new insights rele- 
vant to searches for the elusive CDM particles. This results from the 
unprecedented resolution and convergence (in a dynamical sense) 
of our simulations, as well as from the fact that they provide a sam- 
ple of six Milky Way-like dark matter halos. 

We have measured the probability distribution function of the 
DM mass density between 6 and 12 kpc from the centre of the 
halo, finding it to be made up of two components: a truly smooth 
distribution which scatters around the mean on ellipsoidal shells by 
less than 5% in all the halos of our sample, and a high-density tail 
associated with subhalos. The smooth DM component dominates 
the local DM distribution. With 99.9% confidence we can say that 
the Sun lies in a region where the density departs from the mean 
on ellipsoidal shells by less then 15%. Experimentalists can safely 
adopt smooth models to estimate the DM density near the Sun. 

We find that the local velocity distribution is also expected to 
be very smooth, with no sign of massive streams or subhalo contri- 
butions. The standard assumption of a Maxwellian velocity distri- 
bution is not correct for our halos, because the velocity distribution 
is clearly anisotropic. The velocity ellipsoid at each point aligns 
very well with the shape of the halo. A better fit to the simulations 
is given by a multivariate Gaussian. Even this description does not 
reproduce the exact shape of the distributions perfectly. In partic- 
ular the modulus of the velocity vector shows marked deviations 
from such model predictions. Velocity distributions in our six dif- 
ferent halos share common features with respect to the multivari- 
ate Gaussian model: the low-velocity region is more populated in 
the simulation; the peak of the simulation distribution is depressed 
compared to the Gaussian; at high velocities there is typically an 
excess in the simulation distribution compared to the best-fit multi- 
variate Gaussian. Furthermore the velocity distribution shows fea- 
tures which are stable in time, are reproduced from place to place 
within a given halo, but differ between different halos. These are 
related to the formation history of each individual halo. 

The imprints in the modulus of the velocity vector reflect fea- 
tures in the energy distribution. We explicitly show that the phase- 
space distribution function as a function of energy contains char- 
acteristic wiggles. The amplitude of these wiggles with respect to 
the average distribution function of our sample of six halos rises 
from high to low-binding energies. After appropriate scaling, the 
most bound part of the distribution function looks very similar in 
all halos, suggesting a (nearly) universal shape. The weakly bound 



part of the distribution, on the other hand, can deviate in any given 
halo by an order of magnitude from the mean. 

We have used our simulations to predict detector signals for 
WIMP and axion searches. We find that WIMP recoil spectra can 
deviate about 10% from the recoil rate expected from the best-fit 
multivariate Gaussian model. The energy dependence of these de- 
viations looks similar in all six halos; especially at higher binding 
energies. We find that the annual modulation signal peaks around 
the same day as expected from a multivariate Gaussian model 
with no clear trend over our halo sample for varying recoil ve- 
locity thresholds. The maximum recoil modulation amplitude, on 
the other hand, shows a clear threshold-dependent difference be- 
tween the signal expected for a multivariate Gaussian model and 
that estimated from the simulation. We have also explored the ex- 
pected signal for direct detection of axions. We find the axion spec- 
tra to be smooth without any sign of massive streams. The spectra 
show characteristic deviations from those predicted by a multivari- 
ate Gaussian model; the power at low and high frequencies is higher 
than expected. The most pronounced effect is that the spectra peak 
at lower frequencies than predicted. Since the frequencies in the 
axion detector are directly proportional to the kinetic energy of the 
axion particles, the bumps in the DM velocity and energy distribu- 
tions are also clearly visible in the axion spectra. All the effects on 
the various detector signals are driven by differences in the veloc- 
ity distribution. Individual subhalos or streams do not influence the 
detector signals however, since they are sub-dominant by a large 
factor in all six halos. 

Our study shows that, once direct dark matter detection has 
become routine, the characterisation of the DM energy distribution 
will provide unique insights into the assembly history of the Milky 
Way halo. In the next decade, a new field may emerge, that of "dark 
matter astronomy". 
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